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Differential-Geometric  Techniques  for  Solving  Differential 
Algebraic  Equations 

by 

Florian  A.  Potra1  and  Werner  C.  Rheinboldt2 


1.  Introduction 

Differential-algebraic  systems  of  equations  (DAEs)  arise  in  many  areas  of 
science  and  engineering.  In  particular,  the  equations  of  motion  for  a 
constrained  mechanical  system  considered  in  this  volume  are  usually 
modelled  as  a  second  order  DAE.  In  recent  years  the  literature  on  the 
numerical  solution  of  such  systems  has  grown  rapidly,  (see-e^g.  the  -~>- 
monographs  ft],  [4]).  However,  up  to  now,  existence  theories  for  nonlinear 
DAEs  are  available  only  for  a  few  selected  classes  of  systems. 


In  [10]  a  differential-geometric  approach  was  introduced  for  the  analysis  of 
the  solution  properties  of  a  class  of  linear  DAEs,  and  in  [12]  this  approach  was 
extended  to  general,  semi-implicit,  nonlinear  equations  of  first  and  second 
order;  that  is,  to  systems  with  separated  algebraic  and  differential  equations.  ~  A 

Moreover,  it  was  shown  there  that  these  results  lead  to  a  general  local 
parametrization  approach  for  the  computational  solution  of  these  systems. 

The  differential-geometric  approach  is  based  on  the  observation  that,  as 
long  as  we  expect  the  solution  of  a  DAE  to  be  some  smooth  path  in  the  space 
of  dependent  variables,  the  system  must  define  a  dynamical  system  in 
suitable  subsets  of  that  space.  While  this  connection  with  dynamical  systems  is 
obvious  for  ordinary  differential  equations  (ODEs),  the  same  is  certainly  not 
true  for  DAEs.  In  fact,  besides  [10],  [12]  we  are  only  aware  of  [9]  as  the  only 

1  Department  of  Mathematics,  The  University  of  Iowa,  Iowa  City,  Iowa  52242 
2Dept.  of  Mathematics  and  Statistics,  University  of  Pittsburgh,  Pittsburgh,  PA 
15260.  The  work  of  this  author  was  in  part  supported  under  ONR-grant  N-00014-90- 
J-1025  and  NSF-grant  CCR-8907654. 
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further  study  that  specifically  addresses  this  connection. 

This  article  is  intended  to  be  an  introductory  overview  of  certain  of  the  cited 
differential-geometric  results  and  of  some  of  the  numerical  algorithms  that  can 
be  developed  from  them.  In  order  to  make  the  presentation  widely  accessible, 
we  begin  in  section  2  with  a  brief  summary  of  concepts  and  notations  about 
sub-manifolds  of  finite-dimensional  spaces  and  their  first  and  second  tangent 
bundles.  Then  section  3  presents  an  overview  of  basic  results  about  dynamical 
systems  in  the  form  needed  for  our  DAE-theory.  This  is  followed  in  sections  4 
and  5  by  a  presentation  of  the  principal  results  of  [12]  on  the  existence  and 
uniqueness  of  solutions  and  the  use  of  local  parametrizations.  However,  in 
order  to  avoid  some  of  the  technical  aspects  of  the  general  theory  we  restrict 
ourselves,  for  the  most  part,  to  DAEs  with  certain  linearity  properties  including, 
in  particular,  the  equations  of  motion  for  constrained  mechanical  system 
mentioned  earlier.  In  section  6  we  discuss  the  computational  implementation 
of  local  parametrizations  and  identify  two  useful  parametrizations  induced  by 
the  tangent  space  of  the  constraint  manifold  and  by  a  so-called  "coordinate 
partitioning",  respectively.  Then  section  7  presents  algorithms  based  on  the 
application  of  explicit  and  implicit  multi-step  methods  to  the  resulting  local 
dynamical  systems  and  in  section  8  versions  of  these  algorithms  for  the 
solution  of  the  Euler-Lagrange  equations  are  formulated.  These  algorithms 
generalize  those  previously  discussed  in  [5]  and  [8].  Finally,  in  section  9  we 
give  some  computational  results  for  a  four-bar-linkage  mechanism.  . 


2.  Differential-Geometric  Background 

In  this  section  we  collect  some  basic  material  needed  throughout  the 
remainder  of  the  presentation.  For  further  details  we  refer  to  standard  text  on 
differential  geometry  such  as  [6]  or  [13]. 

We  begin  with  some  standard  terminology.  If  U  is  an  open  set  of  Rn  then  a 

mapping  F:U  ->  Rm  from  U  into  Rm  is  of  class  Cp  ,  p>0,  on  U  if  all  partial 
derivatives  of  F  up  to  and  including  order  p  exist  and  are  continuous  in  U. 
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More  generally,  on  an  arbitrary  set  S  of  Rn  a  map  F:S  ->  Rm  is  of  class  Cp  if  for 
each  x  e  S  there  exists  an  open  set  UcRn  containing  x  and  a  mapping 

G:U  -*  Rm  of  class  Cp  that  coincides  with  F  throughout  UnS.  A  map  F:S  ->  T 
between  the  sets  S  c  Rn  and  T  c  Rm  is  a  homeomorphism  if  F  is  a  one-to-one 
mapping  from  S  onto  T  and  both  F  and  its  inverse  FTT  ->  S  are  continuous.  A 

homeomorphism  F:S  -»  T  is  a  Cp-diffeomorphism  between  S  and  T  and  both  F 
and  F-1  are  of  class  Cp. 

As  usual,  L(Rn,Rm)  is  the  space  of  all  linear  mappings  with  domain  Rn  and 
range  in  Rm  and  by  L2(Rn,Rm)  the  space  of  all  bilinear  maps  from  Rn  to  Rm.  For 
any  Carnap  F:U  c  Rn  Rm,  the  derivative  of  F  at  x  e  U  is  the  linear  map  DF(x) 
in  L(Rn,Rm)  defined  by  DF(x)h  =  lim^g  (F(x+th)  -  F(x)].  In  other  words,  DF(x)  is 
the  linear  map  that  corresponds  to  the  mxn  matrix  of  first  partial  derivatives  of  F 

at  x.  Analogously,  if  F  is  of  class  C2,  then  the  second  derivative  of  F  at  x  is  the 
bilinear  map  D2F(x)  e  L2(Rn,Rm)  defined  by  the  second  partial  derivatives  of  F 
at  x.  Note  that  when  F  is  a  C^diffeomorphism  between  the  open  sets  UcRn 
and  V  c  Rm  then  n=m  and  DF(x)  is  nonsingular  at  all  points  x  e  U. 

A  subset  M  c  Rn  is  a  d-dimensional  Cp-sub-manifold  of  Rn  if  for  each  point 
x  e  M  there  exists  an  open  set  UcRn  containing  x  such  that  the  neighborhood 

U  n  M  of  x  on  M  is  Cp-diffeomorphic  to  an  open  subset  V  of  Rd.  Any  particular 
diffeomorphism  ":U  n  M  ->  V  is  called  a  chart  onllnM  and  its  inverse  a  local 
coordinate  system  on  U  n  M. 

Note  that  by  this  definition  any  open  subset  U  of  Rn  is  an  n-dimensional 
C°°-sub-manifold  of  Rn.  The  tangent  space  TXU  of  this  manifold  U  at  any  point 
x  e  U  is  the  n-dimensional  linear  space  {x}xRn  and  the  first  and  second  tangent 
bundles  TU  and  T^U  of  U  are  the  2n-  dimensions-  and  4n-dimensional 
submanifolds  UxRn  of  (Rn)2  and  Ux(Rn)3  of  (Rn)4,  respectively3. 

The  classical  example  of  a  2-dimensional  C“-sub-manifold  of  R3  is,  of 
course,  the  unit  sphere  S2  =  { (£/n,Q  e  R3;  Z,2  +  ri2  +  X2  =  1  }.  More  generally, 


3We  use  here  the  notation  (Rn)k  =  RnxRn  ...  xRn  (k-times). 
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let  F:E  c  Rn  ->  Rm,  n  >m,  be  some  mapping  of  class  Cp  ,p>1 ,  on  the  open  set  E 
of  Rn.  A  point  x  e  E  is  a  regular  point  of  F  if  dim  DF(x)Rn  =  m;  that  is,  if  the 
derivative  DF(x)  has  full  rank  m.  A  point  be  Rm  is  a  regular  value  of  F  if  the 
inverse  image  F^b)  =  {  x  e  E,  F(x)  =  b}  consists  only  of  regular  points.  Then 
the  following  result  holds: 

Theorem  1 :  Let  F:E  c  Rn  -»  Rm  ,  n  >  m,  be  a  mapping  of  class  Cp  ,p>1 ,  on  the 
open  set  E  of  Rn.  Then  for  any  regular  value  be  Rm  the  inverse  image  F(‘^(b)  is 

either  empty  or  an  (n-m)-dimensional  Cp-sub-manifold  of  Rn. 

From  now  on,  let  F:E  c  Rn  -»  Rm  ,  d=n-m  >  0,  be  a  given  mapping  of  class 

Cp  ,p>3,  on  the  non-empty  open  set  E  of  Rn  such  that  DF(x)Rn  =  Rm  for  xeE. 
Then  each  member  of  the  family  of  sets 

(2.1)  Mb  =  {  x  e  E,  F(x)  =  b},  be  F(E), 

is  a  non-  empty  d-dimensional  Cp-sub-manifold  of  Rn.  Clearly,  any  point  xqsE 
belongs  to  exactly  one  of  these  manifolds,  namely  that  for  b  =  F(x0);  we  shall 
call  this  the  manifold  of  the  family  through  the  point  x0. 

For  any  point  x  e  Mb  the  tangent  space  TxMb  of  Mb  at  x  is  defined  by 

(2.2)  TxMb  =  { (x,p)  e  TxRn;  DF(x)p  =  0  }. 

Clearly,  TxMb  is  a  d-dimensional  linear  subspace  of  the  n-dimensional  linear 
space  TxRn  =  (x}xRn.  The  tangent  bundle  TMb  of  Mb  is  the  disjoint  union  of  all 
tangent  spaces  TxMb  for  x  e  Mb;  that  is, 

(2.3)  TMb  =  { (x,p)  e  TE;  F(x)  =  b,  DF(x)p  =  0  } . 

In  other  words,  the  points  (x,p)  e  TMb  are  the  zeroes  of  the  mapping 


H:  Ex  Rn  ->  Rm  x  Rm;  H(x,p)  =  (  ^*jx'  £  j,  V  (x,p)  e  Ex  Rn 
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Since  DF(x)  was  assumed  to  have  full  rank  on  E  it  follows  that  the  derivative 

DH(,p)=(  0  ) 

l  D2F(x)p  DF(x)  / 

has  full  rank  on  TE  and  hence,  by  Theorem  1 ,  that  the  tangent  bundle  is  a 
non-empty,  2d-dimensional  Cp_1-sub-manifold  of  TRn.  Therefore,  TMb  itself  has 

a  tangent  bundle,  namely,  the  4d-dimensional  Cp-2-sub-manifold  T(TMb)  = 
T2Mbof  the  4n-dimensional  product  space  T2Rn  defined  by 

(2.4)  T2Mb  =  {((x,y),(p,q))  eT2E;  F(x)=b,  DF(x)p  =  0,  DF(x)q+D2F(x)(y,p)=0  }. 


3.  Vectorfields 

A  C°-vectorfield,  o  >  1 ,  on  some  open  subset  E0  of  Rn  is  a  Ca-mapping  on 
E0  such  that 

(3.1 )  k:  E0  ->  TE0  ;  n(x)  =  (x,  9(x)),  V  x  e  E0  . 

An  integral  curve  of  n  through  a  point  x0  e  E0  is  any  C°-path  E,:  J  ->  E0,  defined 
on  some  open  interval  J  <=  R1  containing  0  e  J,  for  which  (£(t),£'(t))  =  n(E,(t))  for 
t  e  J  and  E;(0)  =  x0  ;  that  is,  which  solves  the  initial  value  problem 

(3.2)  x'  =  8(x),  x  e  E0  ,  x(0)  =  x0. 

As  before,  let  F:E  cR"^  Rm,  d=n-m  >  0,  be  a  given  mapping  of  class  Cp, 
p>3,  on  the  non-empty  open  set  E  of  Rn  such  that  DF(x)Rn  =  Rm  for  x  e  E,  and 

consider  the  family  (2.1)  of  d-dimensional  Cp-sub-manifolds  of  Rn.  We  call  the 
vectorfield  (3.1)  tangential  to  this  fami'y  of  manifolds  (2.1)  if  rc(x)eTxMb  for  all 
xe  E0nMb.  Since  xe  Mb,  b=F(x),  for  any  xe  E0nE,  this  requires  that 


(3.3) 


DF(x)0(x)  =0,  V  x  e  E0nE. 
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Thus,  for  any  integral  curve  %:  J  ->  E0nE  of  k  through  x0  e  E0nE  it  follows  that 
DF(c,(t))£,,(t)  =  0,  t  e  J.  Since  £(J)  is  a  connected  subset  of  the  open  set  E0nE, 
the  integral  mean-value  theorem  guarantees  that 

F($(t))  -  F(^(0))  =  f  DF(£(s))^(s)ds  =  0,  Vte  J. 

Jo 

In  other  words,  the  path  J  E0nE  remains  on  the  manifold  Mb,  b  =  F(x0), 
through  the  initial  point  x0,  and  (3.3)  implies  that  (£(t),  £'(t))  e  TMb,  fortej. 

From  this  and  the  standard  theory  of  ODEs  (see  e.g.  [6]  or  [13])  we  obtain 
now  the  following  existence  and  uniqueness  result: 

Theorem  2:  Let  (3.1 )  be  a  C°-vectorfie!d,  a  > 1 ,  on  the  (non-empty)  open 
subset  E0  of  Rn.  Then  the  following  results  hold: 

(i)  There  exists  a  C°-integral  curve  J  ->  E0  of  it  through  each  x  e  E0 
defined  on  an  open  interval  J  containing  0.  Moreover,  any  two  such  curves  are 
equal  on  the  intersection  of  their  domains. 

(ii)  If  k  is  tangential  to  the  family  (2.1)  and  E0nE  is  non-empty,  then  any 
integral  curve  £,:J  ->  E0r»E  of  rc  through  x0  e  E0nE  satisfies  (E,(t),  ^'(t))  e  TMb, 
b=F(x0),  for  all  t  e  J. 

(iii)  The  union  of  the  domains  of  all  integral  curves  of  k  through  a  point 

xe  U0  is  an  open,  possibly  unbounded  interval  Jx  =  (x.(x),  x+(x)).  There  exists  a 

C°-integral  curve  £*:JX  ->  E0,  of  n  through  x,  and  Jx  is  the  largest  interval  on 
which  such  an  intergal  curve  exists. 

(iv)  If  x+(x)  <  for  some  x  e  E0,  then  for  any  compact  set  CcE0  there  exists 
a  5  >  0  such  that  ^*(t)  €  C  for  t  >  x+(x)-5.  A  corresponding  result  holds  when 
x(x)  >  -oo. 

(v)  The  set  D{n)  =  {(t,x)  e  R^Uq;  te  Jx}  is  open  in  R’xEq  and  contains 
{0}xU0.  Moreover,  the  global  flow  y.D( n)  -»  E0  ,Y(t,x)=^*(t),  te  Jx ,  of  k  is  of 
class  C°  on  D{n). 


Consider  now  a  second  order  initial  value  problem  x"  =  0(x,x'),  x  e  E0, 
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x(0)=x0,  x'(0)=y0,  which,  of  course,  may  be  written  in  the  first  order  form 

(3.4)  x'  =  y,  y'  =  0(x,y),  x  e  E0,  x(0)  =  x0,  x’(0)  =  y0. 

Thus,  we  encounter  in  this  case  a  vector-field  of  the  form 

(3.5)  n:  En  =  E0xY cTE0-»T2E0  ;  *(x,y)  =  ((x,y),(y,e(x,y))),  V  (x,y)  e  E0xY 

on  a  subset  of  the  tangent  bundle  TE0  of  E0.  Note  that  the  second  and  third 
component  of  the  image  vector  are  identical.  In  other  words,  (3.5)  represents  a 
sub-class  of  the  vector  fields  on  En  called  the  vector  fields  on  En  that  are 
consistent  with  a  second  order  ODE.  We  assume  always  that  the  domain  En  of 
n  is  some  open  subset  of  TE0. 


An  integral  curve  of  (3.5)  through  a  point  (x0,  y0)  e  En  is  now  a  C°-path 
^:J-»  E0,  defined  on  some  open  interval  JcR1  containing  0  e  J,  for  which 


(3.6) 


memjeE, ,  ^oh^,  ?(<>)-y0 
(mmumzm = nmsm  v  t  e  j; 


that  is,  which  is  a  solution  of  the  initial  value  problem  (3.4). 


As  before,  the  vectorfield  (3.5)  is  called  tangential  to  the  family  of  tangent 
bundles  TMb,  be  F(E),  of  (2.1 )  if  ji(x,y)  eT(x  y)(TMb)  for  (x,y)e  En  nTMb;  that  is,  if 

(3.7)  DF(x)0(x,y)  +  D2F(x)(y,y)  =  0,  V  (x,y)  e  En  n(ExkerDF(x)) 

Then,  for  any  integral  curve  ^:J  E0nE,  of  n  through  (x0,y0)  e  E^nTMh  such 
that  DF(^(t))^’(t)  =  0,  t  e  J,  the  integral  mean-value  theorem  provides  that 


DFfc(t))?(t)-DF(x0)yo 


IDF(5(s))4"(s)+D2F(%(s))(£,'(s).£.-(s))3ds  =  0.  V  t  e  J 

0 
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Hence,  the  path  £:J  ->  E0nE  remains  on  the  tangent  bundle  TMb>  b  =  F(x0), 
corresponding  to  x0  and  (3.7)  implies  that  ((^(t),^,(t)),(^’(t),^"(t)))  eT2Mb,  fort  e  J. 

Again  an  existence  and  uniqueness  result  of  the  form  of  Theorem  2  holds. 
We  note  here  only  the  following  shortened  version. 

Theorem  3:  Let  (3.5)  be  a  C°-vectorfield,  o  >  2,  on  the  (non-empty)  open 
subset  En  <z  TE0  cz  TRn.  Then  there  exists  an  integral  curve  J  U0  of  n 
through  each  (x,y)  e  En  defined  on  an  open  interval  J  containing  0.  Moreover, 
any  two  such  curves  are  equal  on  the  intersection  of  their  domains.  If  k  is 
tangential  to  the  family  (2.1)  and  E0nE  is  non-empty,  then  any  integral  curve 
^:J  -»  E0nE  of  tc  through  (x0,y0)  e  EnnTMb  ,  b=F(x0),  for  which  DF(£(t))£'(t)  =  0 
for  t  in  J,  satisfies  ((^(t),^'(t)),(4'(t)£"(t)))  eT2Mb,  for  all  t  e  J. 

4.  A  Class  of  First  and  Second  Order  DAEs 

In  this  section  we  begin  with  the  autonomous,  first  order  DAE 

(4.1)  Ft(x)  =  b 

F2(x,  x',  z)  =0 
for  which  we  assume  that 

(4.2a)  F, :  ExcRn  -»  Rr  and  F2:  E2  =  ExxEpxEz  c  (Rn)2  xRm  Rs  are 

of  class  CP,  p  £  3,  on  their  domains,  where  ExcRn,  EpcRn,  and 
Ez  c  Rm,  are  non-empty  open  sets,  and  r<n<r+s=n+m;  and 
(4.2b)  for  each  (x.p.z)  e  E2  the  matrix 

I  DFt(x)  0  \ 

|  DpF2(x,p,z)  DzF2(x,p,z)  J 

is  non-singular. 

From  (4.2b)  it  follows  that  DF^xJR"  =  Rr  for  all  x  e  Ex  and  hence,  as  we  saw 
in  section  2,  that  each  member  of  the  family  of  sets 
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(4.3)  Mb  =  {xe  Ex,F1(x)  =  b},  be  F  1(EX), 

is  a  non-  empty  (n-r)-dimensional  Cp-sub-manifo!d  of  Rn.  For  any  x0  e  Ex  we 
call  the  unique  manifold  (4.3)  with  b  =  F.,(x0)  the  constraint  manifold  through 
that  point. 

For  given  be  F^EJ  a  Ca-solution,  1  <  a  £  p— 1 ,  of  (4.1)  consists  of  two  C°- 
paths  x:  J  ->  Ex  and  z:  J  -4  Ez,  defined  on  some  open  interval  J  of  R1,  such 
that 

(4. .4)  (x(t),x'(t),z(t))  e  E2  ,  F,(x(t))  =  b,  F2(x(t),x'(t),z(t))  =  0,  V  t  e  J. 

This  necessitates  that  DF.,(x(t))x'(t)  =  0,  V  t  e  J.  Hence,  if  (x0,p0,z0)  e  E2  is  any 
point  on  a  solution  (4.4);  that  is,  if  x(t0)  =  x0,  x'(t0)  =  p0,  z(t0)  =  z0  for  some  t0  e  J, 
then  we  must  have  DF^XqJpq  =  0,  and  F2(x0,p0,z0)  =  0,  while,  of  course, 
automatically  x0  e  Mb  for  bsF^Xo). 

This  suggests  the  definition  of  the  initial  data  map 

(4.5)  H:  E2  -»  RrxRs,  H(x,p.z)  =  (  ),  (x,p.z)  s  E2 . 

Evidently,  for  any  given  (x,p,z)  e  E2,  the  partial  derivative  Dp  ZH  of  H  with 
respect  to  p  and  z  is  the  matrix  in  (4.2b).  By  Theorem  1  this  implies  that 

(4.6)  K  =  { (x,p,z)  e  EP  ;  H(x,p,z)  =  0  } 

is  a  n-dimensional  CP_1-submanifold  of  (Rn)2xRm,  the  initial  data  manifold  of 
the  problem. 

A  general  existence  and  uniqueness  theory  for  (4.1)  was  given  in  [12]  and 
requires  a  closer  analysis  of  the  relationship  between  the  initial  data  manifold 
K  and  the  constraint  manifolds  (4.3)  utilizing  the  theory  of  covering  spaces.  For 
the  applications  considered  here  it  suffices  to  restrict  attention  to  the  case 
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when  F2  is  linear  in  p  and  z;  that  is,  when  (4.1)  has  the  special  form 

(4.7)  F1  (x)  =  b 

A(x)x  +  B(x)z  +  G(x)  =  0, 


Here  (4.2a)  requires  that  F1 :  ExcRn->  Rr,  A:Ex->  L(Rn,Rs) ,  B:  Ex->  L(Rn,Rs),  and 

G:  Ex->  Rs  are  of  class  Cp,  p  >  3,  on  the  open  non-empty  set  ExcRn,  and  that 
r<n<r+s=n+m.  Moreover,  condition  (4.2b)  is  equivalent  with  the  assumption 
that  for  each  x  e  Ex  the  matrix 


(4.8) 


I  DFi(x)  0  \ 
1  A(x)  B(x)  I 


is  non-singular. 


Hence,  in  this  case  the  equation  H(x,p,z)  =  0  has  for  each  x  e  Ex  a  unique 
solution 

,49,  fP  0  YY  0  >1 

'  '  '  U  J  U(x)  J  l  A(x)  B(x)  J  (-G(x)J 

Then 

(4.10)  tc:  Ex  ->TEX,  n(x)  =  (x,  T](x)),  xe  Ex 


defines  a  Cp-vectorfield  on  the  open  set  Ex  which  by  definition  of  the  initial 
data  manifold  K  is  tangential  to  the  family  of  constraint  manifolds  (4.3). 

Now,  consider  (4.7)  together  with  the  initial  condition  x(0)  =  x0  e  Ex  and  set 
p0=r[(x0),  z0=£(x0)  which  implies  that  (x0,  p0,  z0)  e  K.  Then  Theorem  3  holds 
and  there  exists  an  integral  curve  ^:J  -»  Ex  of  n  through  x0  for  which 
(£(t),£’(t))e  TMb,  b=F(x0),  for  all  t  €  J.  In  other  words,  the  paths  x:J-»Ex ,  z:J-»Ez, 
defined  by  x(t)=^(t),  z(t)=C(^(t)),  te  J,  constitute  a  solution  of  (4.7). 

Instead  of  repeating  all  of  Theorem  3  we  summarize  this  result  only  as 
follows: 
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Theorem  5:  If  (4.2a/b)  holds  for  (4.7),  then,  for  any  given  x0  e  Ex ,  there  exists  a 
unique,  maximally  extended  Cr-solution  x:  J  ->  E1  and  z:  J  — >  Ez,  of  (4.7)  which 
satisfies  the  initial  conditions  x(0)  =  x0,  x'(0)  =  p0  =  r|(x0),  z(0)  =  z0=^(x0). 

As  noted,  a  corresponding  result  for  the  general  equation  (4.1)  was  proved 
in  [12].  There  the  theory  was  also  extended  to  the  second  order  system 

(4.11)  F^x)  =  b 

F2(x,  x',  x",  z)  =  0 

subject  to  the  conditions 

(4. 1 2a)  F1 :  ExcRn  ->  FT  and  F2:E2  =  ExxEyxEqxEz  c  (Rn)3  xRm  -►  Rs 

are  of  class  CP,  p  >  4,  on  their  domains  where  Ey,Eu,E;,cRn, 

x  y  z 

and  E2  c  Rm,  are  non-empty  open  sets  and  r<n<r+s=n+m; 
(4.12b)  for  each  (x,y,q,z)  g  E2  the  matrix 

I  DF^x)  0  \ 

lDqF2(x,y,q,z)  DzF2(x,y,q,z)| 

is  non-singular. 

As  before  we  see  that  each  member  of  the  family  of  sets  (4.3)  is  a  non¬ 
empty  (n-r)-dimensional  Cp-sub-manifold  of  Rn,  the  constraint  manifold  of 
(4.1 1 )  through  x0. 

It  is  natural  to  reduce  (4.11)  to  the  first  order  system: 

F,(x)  =  b 

(4.13)  x’-y  =0 

F2(x,  y  v’  z)  =  0. 

With  the  combination  (x,y)  as  new  differential  variable,  this  constitutes  a  DAE  of 
the  form  (4.1 )  for  which  (4.2a)  turns  out  to  be  valid.  However,  (4.2b)  does  not 
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hold  since  the  corresponding  matrix 

'  DF-i(x)  0  0 

In  0  0 

0  DqF2(x,y,q,2)DzF2(x1y>q,z)/ 

is  singular.  This  failure  of  (4.2b)  is  hardly  surprising  since  we  should  expect 
(4.1 1 )  to  induce  a  vector  field  on  the  tangent  bundle  TEX  that  is  consistent  with 
a  second  order  equation. 

For  any  b  e  Fi(Ex)  a  Ca-solution,  1<c<p-2,  of  (4.1)  consists  of  C°-  paths 
x:  J  Ex  and  z:J  — >  E2,  defined  on  some  open  interval  J  of  R1,  such  that 

(4.14)  (x(t),x’(t),x"(t),z(t))  e  E2,F1(x(t))  =  b,  F2(x(t),x’(t),x"(t),z(t))  =  0,  VteJ. 

This  implies  that  DF1(x(t))x'(t)  =  0,  DF^x^x’^t)  +D2F1(x;(x’(t),x'(t))=0,  for  all  t 
in  J,  which  means  that  y:J->T2Mb  ,  y(t)  =  ((x(t),x'(t)),(x'(t),x"(t)),  t  e  J,  is  a  path 
on  T2Mb.  In  analogy  to  the  first  order  case,  this  suggests  the  definition 

(4.15)  H:E2  c  (Rn)3xRm  -»  Rr+S 

H(x,y,q,z)  =  (DF1(x)q+D2F1(x)(y,y),  F2(x,y,q,z))  ,  V  (x,y,q,z)eE2 

for  the  initial  data  map  of  (4.1 1 ).  Once  again,  (4.1 2b)  implies  that  the  solution 
set 

(4. 1 6)  K  =  { (x.y.q.z)  e  E2  ;  H(x,y,q,z)  =  0 } 

is  an  n-dimensional  Cp_2-submanifold  of  (Rn)3xRm,  called  the  initial  data 
manifold  of  (4.1 1). 

As  before,  we  shall  not  present  here  the  general  theory  but  restrict 
ourselves  to  the  linear  case 

F  i  (x)  -b 

A(x,x')x"  +  B(x,x')z  +  G(x,x')  =  0, 


(4.17) 
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Then  (4.12a)  requires  that  the  maps  F, :  ExcRn  Rr ,  A:  ExxEy-»  L(Rn,Rs) , 
B:ExxEy  L(Rm,Rs)  and  G:  ExxEy  Rs  are  of  class  CP,  p>4,  on  their  domains 
where  Ex  ,Eyc  Rn  are  non-empty  open  sets  and  r<n^r+s=n+m.  Moreover, 
(4.12b)  is  equivalent  with  the  assumption  that  the  matrix 


(4.18) 


(  DF i  (x)  0  ) 

1  A(x,y)  B(x,y)  | 


is  non-singular  for  all  (x,y)  e  ExxEy. 


Then,  as  in  the  first  order  case,  the  equation  H(x,y,q,z)  =  0  has  for  each  (x.y) 
in  ExxEy  a  unique  solution 


(4.19) 


'  n(x,y)  ' 

)  =  Jf  DFi(x)  0  ^ 

-i  f  D2F,(x)(y,y) ' 

J  l 

,  C(x,y)  , 

j  l  A(x,y)  B(x,y)  J 

,  G(x,y)  , 

and  hence 


(4.20)  k:  En  =  ExxEy  c  TEX  T^E, ,  rc(x)  =  ((x,y),(y,Ti(x,y))),  (x,y)  e  En 

defines  a  Cp_1-vectorfield  on  the  open  set  En  of  the  tangent  bundle  TEX. 
Clearly,  by  definition  of  the  initial  data  manifold  K,  this  vectorfield  is  tangential 
to  the  family  of  tangent  bundles  TMb,  be  F^E^  of  the  constraint  manifolds  (4.3) 
of  (4.1 7);  that  is, 

(4.21)  DF-i  (x)rj(x,y)  +  D2F1(x)(y,y)  =  0,  V  (x,y)  e  Exx(Eyn  kerDF^x)) 
Now,  consider  (4.17)  together  with  the  initial  condition 

(4.22)  x(0)  =  x0  ,  y(0)  =  y0  ,  (x0,y0)  e  Ex  x  (Eyn  kerDF, (x0)) 

and  set  q0=Tl(x0,y0),  z0=C(x0,y0)  which  implies  that  (x0,y0,  q0,  z0)  e  K.  Thus  by 
Theorem  3  there  exists  an  integral  curve  £:J  -»  Ex  of  it  for  which  (3.6)  holds. 
By  definition  of  n  we  have  ((^(t),^,(t)),(^,(t),^"(t)))  e  K,  t  e  J,  and  hence,  in 
particular  DF^ft))?^)  =  0,  t  e  J.  Thus  it  follows  that  ((£(t)£'(t)).(£'(t),£"(t)))  e 
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T2Mb,  b=F.,(x0),  and  therefore  that  the  paths  x:J-»Ex  ,  z:J->Ez,  defined  by 
x(t)=£(t),  z(t)=C(^(t)),  te  J,  constitute  a  solution  of  (4.7). 

Once  again,  instead  of  repeating  all  of  Theorem  4  we  summarize  this  result 
only  in  the  brief  form: 

Theorem  6:  If  (4.1 2a/b)  holds  for  (4.1 7),  then,  for  any  given  initial  data  (4.22), 
there  exists  a  unique,  maximally  extended  Cp-solution  x:  J  -»  Ex ,  z:  J  Ez,  of 
(4.17)  which  satisfies  the  initial  conditions  x(0)  =  x0,  x'(0)  =  y0,  x"(0)=ri(x0,y0), 
z(0)  =  C(x0,y0). 

For  a  corresponding  result  for  the  general  equations  (4.1 1 )  see  again  [1 2]. 


5.  Local  Parametrizations 

As  shown  in  [12]  the  results  about  the  semi-implicit  DAEs  (4.1)  and  (4.11) 
lead  to  a  general  local  parametrization  approach  for  the  numerical  solution  of 
these  systems.  Once  again,  we  shall  restrict  our  discussion  here  to  the  linear 
systems  (4.7)  and  (4.17). 

Consider  first  the  system  (4.7),  and  suppose  that  we  are  in  the  setting  of 
Theorem  5.  Thus,  we  wish  to  compute  a  numerical  approximation  of  the 
unique  Cp-solution  x:  J  -»  Ex ,  z:  J  ->  Ez  of  (4.7)  which  satisfies  the  initial 
conditions  x(0)  =  x0,  p0  =  n(x0),  z0  =  £(x0)  for  given  x0  e  Ex.  Here  q:Ex  ->  Ep  and 
£:EX  Ez  are  the  mappings  (4.9).  As  we  saw,  the  problem  of  solving  (4.7)  in 
Ex  then  reduces  to  that  of  solving  the  explicit  initial  value  problem 

(5.1)  x’=  r|(x),  x  e  Ex,  x(0)  =  x0. 

Since  the  desired  solution  x:  J  ->  E0  of  (5.1)  through  x0  remains  on  the  tangent 
bundle  TMb  of  the  constraint  manifold  Mb  ,  b  =F(x0),  through  x0  it  is  natural  to 
work  with  a  local  coordinate  system  on  Mb. 
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Consider  any  point  xc  e  Mb,  on  the  manifold  Mb  through  x0.  The  subscript  'o' 
indicates  here  that  xc  will  later  represent  a  'current'  (approximate)  point  on  the 
solution.  Now  choose  a  linear  subspace  T  c  R",  dim  T  =  n-r,  such  that  T1n  ker 
D^Xg)  =  {0}.  Then  it  follows  from  the  implicit  function  theorem  that  there  exist 
open  neighborhoods  V  of  the  origin  in  T  and  Sx  of  xc  in  Rn,  respectively,  as 

well  as  a  Cp-map  wiV-VFk:  Rn  with  w(0)  =  0,  such  that  MbnSx=vF(V)  where  *¥ 
is  the  local  coordinate  map 

(5.2)  *F:  V  ->  Rn,  'F(t)  =  xc  + 1  +  w(t),  t  e  V. 

In  practice  it  is  often  useful  to  work  with  local  coordinate  mappings  that 
correspond  to  the  tangent  spaces  of  Mb;  that  is,  to  choose  T=kerDF(xc).  Note 
that  then  Dw(0)  =  0,  (see,  e.g.,  [11]). 

For  the  computation  we  introduce  orthonomal  bases 

Q1  e  L(Rn  r,Rn),  Q2  e  L(Rr,Rn), 

(5.3)  Q1  Rn  r  =  T ,  Q2Rr  =  T1 

Q1*Q1  =  ln,,Q2*Q2=lr  ,Q2*Q1  =  0 

for  T  and  T1.  Then  we  have  the  transformed  local  coordinate  map 

(5.4)  <J>:  U  =  Cyv  ->  Rn,  O(u)  =  xc  +  C^u  +  Q2co(u),  u  e  Uc  Rn'r 

where  to(u)  =  (Q2|T-L)'1w((Q1|T)'1  u),  u  eU,  is  again  of  class  Cp  and  satisfies 
co(0)=0.  From  Ft(0(u))  =  0  we  find  that  DO(u)Rn'r  c  ker  DF1(0(u))  for  u  e  U 
and,  by  (5.3),  that 

(5.5)  Q1*DO(u)  =  Q1*[Q1+Q2Dw(u)]  =  ln.r 
whence,  by  a  dimensionality  argument, 

(5.6)  DO(u)Rn'r  =  ker  DF.,(<I>(u)),  V  u  e  U. 

This  implies  that  for  any  u  e  U  and  c  e  ker  DF1(0(u))  the  equation  DO(u)a=c 
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has  a  solution  a  e  Rn*r  for  which  a  =  Qt*c  and  hence  which  is  unique. 

In  our  local  coordinate  system  we  obtain  for  the  ODE  (5.1)  the  local 
representation 

(5.7a)  D$(u)u'  =  ti(0(u)),U6  U 

which,  because  of  r|(<t>(u))  e  ker  DF^Ofu)),  can  be  written  in  the  form  of  the 
(n-r)-dimensional  explicit  system 

(5.7b)  u’  =  Q,*  ti(<D(u)),  u  e  U. 


Obviously,  with  p  =  Q1ri1  +  Q2Tl2,  e  Rn‘r',  rl2  e  Rr* the  equation  H(x,p,z)  =  0 
becomes  the  non-singular,  linear  system 


(5.8) 


DF^xJQ!  DFi(x)Q2  0 
A(x)Qi  A(x)Q2  B(x) 


Hence,  if,  for  given  x  =  <D(u),  ue  U,  the  solution  of  (5.8)  has  been  found  then  we 
have  Q/  ri(0(u))  =  t\^  and  ^(O(u))  =  z. 


This  local  parametrization  can  also  be  carried  over  to  the  second  order 
system  (4.17).  As  before,  suppose  that  we  are  in  the  setting  of  Theorem  6. 
Thus,  we  wish  to  compute  a  numerical  approximation  of  the  unique  CP'1- 
solution  x:  J  -»  Ex ,  z:  J  ->  Ez  of  (4.17)  which  satisfies  the  initial  conditions 
x(0)=xo,  x'(0)=y0,  x"(0)=r)(x0,yo),  z(0)=C(x0,y0)  where  (x0,y0)e  ExxEy  n  TMb  , 
b=F1(x0),  are  given  and  q:  En  =  ExxEy  Eq  ,  ExxEy  ->  Ez  ,  are  the  mappings 
(4.19).  As  we  saw,  the  problem  of  solving  (4.17)  reduces  to  that  of  solving  the 
explicit  system 


(5.9)  x'=  y,  y’  =  Ti(x,y),  (x,y)  e  En ,  x(0)=x0,  y(0)=y0, 


Since  the  desired  solution  x:  J  ->  Ex  of  (5.1)  through  x0  remains  on  the  second 
tangent  bundle  T^Mb  of  the  constraint  manifold  Mb  through  x0  we  are  led  to 
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working  with  a  local  coordinate  system  on  TMb. 

Once  again  let  xc  be  a  'current'  point  on  the  manifold  Mb  through  x0  and 
choose  a  linear  subspace  T  c  Rn,  dim  T  =  n-r,  such  that  TVikerDF^g)  =  {0}. 
Then,  with  basis  maps  Q1(Q2  satisfying  (5.3),  the  local  coordinate  map  (5.4)  is 
well  defined  and 

(5. 1 0)  @:UxRn‘r  ->  TMb,  0(u,v)  =  (<D(u),D<D(u)v),  u  e  U,  v  e  Rn‘r 

defines  a  local  coordinate  system  on  TMb.  Clearly,  we  can  choose  some 
neighborhood  S0  of  the  origin  of  UxRn'r  such  that  0(u,v)  belongs  to  En  for  all 
(u,v)  in  S0. 

Thus  with  x  =  O(u)  and  y  =  DO(u)v  =  Q1v+Q2Dco(u)v,  the  differential 
equations  (5.9)  assume  the  local  form 

(5.11)  DO(u)u'=  Q.,v+Q2Dgj(u)v 

DO>(u)v'  +  D20(u)(u',v)  =  ri(<t>(u),  C^v+C^DcofuJv). 

Recall  that  for  any  d  e  ker  DF1(0>(u))  the  unique  solution  of  DO(u)a  =  d  is 
a=Q1*d.  Hence,  since  Q^+QgDto^Jv  e  kerDFi(x),  the  first  equation  reduces 
simply  to  u'=  v  and,  with  this,  the  second  equation  can  be  written  as 

D<t>(u)v’  =  •n(x.y)  -D2<D(u)(v,v)  ,  x  =0(u),  y  =  Q1v+Q2Dto(u)v. 

From  DFt  (0(u))D0(u)  =  0,  it  follows,  together  with  (4.21 )  and  (5.9),  that 

(5.12)  0  =  DF^xJD^OjHv.v)  +  D2Fi(x)(DO(u)v,DcI>(u)v) 

=  DF^xJD^^J^.v)  +  D2F1  (x)(y,y) 

=  DFt (x)[D2cD(u)(v,v)  -  r|(x,y)]  ,  V  (u,v)e  S0 

and,  therefore  that,  as  desired,  Ti(x,y)-D2<t>(u)(v,v)  e  ker  DF^x). 

With  this  we  have  shown  now  that  on  S0  the  equations  (5.1 1 )  can  be 
reduced  to  the  explicit  form 
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(5.13)  u'  =  v 

v'  =  Q/  [ti(<1>(u),Q1v+Q2Dco(u)v)  -  D2<J>(u)(v,v)]. 

From  (5.5)  it  follows  that  Q.,*D2<I>(u)=0  ,  so  that  (5.13)  reduces  to 

(5.14)  u'=v 

v'=Q1*ti(0(u),Q1v+Q2Dq)(u)v)  =  Q1*t](<D(u),DO(u)v) 

Now  ,  with  ri(x,y)=Q1ri1  +  Q2Tl2.  e  Rn'r’,  t|2  e  Rr,  the  equation  H(x,y,q,z)  =  0 
becomes  the  non-singular,  linear  system 


f  DFUxJQi 

DF^xJQz  0 

'  D2F,(x)(y,y)  ' 

l  A(x,y)Qi 

A(x,y)Q2  B(x,y) 

AzJ 

,  G(x,y)  , 

Hence,  if,  for  given  x  =  O(u),  y  =  DC>(u)v,  (u,v)  €  S0,  a  solution  of  (5.15)  has 
been  obtained  then  v’  =  Q/ti(<I>(u),D<I>(u)v)  =  ri1  and  z  =^(0(u),D0(u)v). 


6.  Computational  Implementation  of  the  Local  Parametrization 


In  the  previous  section  we  have  used  the  local  coordinate  map  (5.4)  to 
reduce  the  semi-implicit  DAEs  (4.7)  and  (4.17)  to  the  explicit  systems  of  ODEs 
(5.7b)  and  (5.14),  respectively.  Theoretically  any  ODE  solver  can  be  applied  to 
these  ODE  systems  to  obtain  a  local  solution  of  the  original  DAE.  Before 
discussing  this  in  more  detail,  we  consider  first  how  the  local  coordinate  map 
O  of  (5.4)  at  the  current  point  xc  of  the  constraint  manifold  Mb,  b  =  F^Xq),  may 
be  implemented  and  how  we  might  choose  the  matrices  Q1  and  Q2  that  define 
this  local  parametrization. 


By  definition,  for  any  given  u  e  U,  the  point  x  =  O(u)  on  Mb  is  the  solution 
of  the  augmented  system 


(6.1) 


'  F,(x)  Wb 

.  Oi(x  -  Xc)  J  u 


Here,  of  course,  the  matrices  Q1(  Q2  are  assumed  to  satisfy  (5.3)  and  the 
subspace  T=Q1*Rrvr  has  to  induce  a  local  parametrization  of  Mb  at  xc;  that  is, 
we  have 

(6.2)  (Q/  Rn'r  )xn  ker  DF^xJ  =  {0}. 

It  is  readily  seen  that  (6.2)  holds  if  and  only  if  the  Jacobian 


DFt(x) 

Qi 


of  the  mapping  on  the  left  side  of  (6.1)  is  nonsingular  at  xc  which,  in  turn  is 
equivalent  with  the  nonsingularity  of  DF1(xc)Q2. 

As  proved,  for  instance,  in  [7],  under  this  nonsingularity  assumption  there 
exists  an  open  set  U^U  containing  the  origin,  such  that  for  any  u  in 
Newton's  method  applied  to  (6.1)  and  started  from  y<°>  =  xc  +  C^u  produces  a 
sequence  of  points  {y(k)}  that  converges  Q-quadratically  to  x=<I>(u).  A  step  of 
this  process  here  has  the  algorithmic  form: 

(1)  Evaluate  a1  =  F,(yW)  -  b,  a2  =  Q1*(y(k)  -  xc); 

(2)  solve  the  linear  system 


DF,(yW) 


(3)  Set  y(k+1>  =  y(R)-  s. 


-(£> 


Observe  that  for  our  choice  y(0)  =  xc  +  QjU  of  the  starting  point  the  iterates  y(k) 
remain  in  the  affine  subspace  xc  +  Q2Rr. 


Let  s1  =  Q/s,  s2  =  Q2*s,  then,  by  (5.3),  we  have  5  =  0^  +  Q2s2  and  we 
obtain  the  solution  of  the  linear  system  in  step  (2)  by  setting  5^2  and  solving 
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the  (n-r)-dimensional  linear  system 

(6.4)  [DF1(yW)Q2]  s2  =  a,  -  DF^xW)  Qia2 . 

Thus,  when  the  Jacobian  of  F.,(x)  is  readily  available,  then  Newton’s  method 
provides  an  efficient  numerical  procedure  for  implementing  the  local 
coordinate  map  (5.4). 

We  turn  now  to  the  construction  of  the  matrices  Q1t  Q2  for  which  (5.3)  and 
(6.2)  are  satisfied  .  Two  choices  appear  to  be  particularly  useful  in  practice. 

The  first  one  was  already  mentioned  in  the  previous  section  and  defines  the 
parametrization  by  using  the  space  T  =  ker  DF^)  corresponding  to  the 
tangent  space  of  Mb  at  xc.  In  this  case  the  matrices  Qv  Q2  may  be  obtained,  for 
instance,  by  performing  a  QR-factorization  of  the  transposed  Jacobian  of  F1 

(6.5)  DF,(xo)‘-C>(£) 

where  Q  is  an  nxn  unitary  matrix  and  R  an  rxr  upper  triangular  matrix  which,  by 
(4.2b)  is  non-singular.  In  fact,  it  follows  easily  that  then  the  local 
parametrization  is  defined  by  the  matrices 

(6.6)  Q2  =  Q(ev  ...,er) ,  =  Q(er+^  , ... ,  en) 

formed  from  the  first  r  and  last  n-r  columns  of  Q,  respectively.  Here  e1 . en 

denote  the  natural  basis  vectors  of  Rn. 

This  "tangent  space"  parametrization  may  be  somewhat  costly  to  use.  To 
reduce  the  cost  it  is  advantageous  to  choose  a  local  coordinate  space  T  which 
is  spanned  by  n-r  suitably  selected  natural  basis  vectors  of  Rn.  If,  say, 

(6.7)  T1  =  span  { e-h . ejf},  T  =  span  { ek), . ejn } 

then  the  condition  Txnker  DF(xc)  =  {0}  requires  that  we  select  the  indices 

j1 . jr  such  that  the  corresponding  columns  of  DF(xc)  are  linearly 

independent.  It  is  well-known  that  such  a  choice  may  be  derived,  for  instance, 
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from  a  Jordan  factorization  with  row  and  column  pivoting,  or  from  a  singular 
value  decomposition  (SVD)  of  DF^xJ  (see.e.g.  [3]).  The  matrices  of  the  local 
parametrization  induced  by  (6.6)  then  are 

(6.8)  Qi  =  ( ejM . ejn ),  Q2  =  (ei, . ejf) 

and  satisfy  (5.3)  by  construction. 

The  choice  of  the  permutation  [j1(...,jn]  in  (6.7)  and  (6.8)  "partitions"  the 
components  of  any  vector  x  into  a  vector  of  "  independent  coordinates 

"  and  a  vector  =  Q2x  of  "dependent  coordinates  ".  Thus  it  is  natural  to  call 
(6.8)  a  local  parametrization  by  "coordinate  partitioning". 


7.  Explicit  and  Implicit  Multistep  Methods, 

Consider  again  the  first-order  DAE  (4.7)  under  the  assumptions  (4.2a/b) 
which,  by  Theorem  5,  guarantee  the  existence  of  a  solution  x:J  -»  Ex ,  z:J  ->  Ez 
that  satisfies  the  initial  conditions  x(0)  =  x0,  x'(0)  =  p0  =  ti(x0),  z(0)  =  z0=£(x0) 
corresponding  to  the  given  point  x0  e  M„  on  the  constraint  manifold. 

Suppose  that  xc  =  x(tc)  e  Mb  is  a  'current'  point  on  this  solution  where  a 
local  parametrization  (5.4)  has  been  constructed.  In  some  region  near  xc  we 
wish  to  compute  a  sequence  of  approximation 

(7.1 )  Xj  =  x(tj),  Zj  =  z(tj),  tj=  tc+  ih  ,  i=1 ,2 . N 

of  points  x(tk)  along  the  solution.  More  precisely,  we  suppose  that  the 
constants  h  and  N  are  chosen  such  that  the  points  x(tk)  belong  to  the 
neighborhood  of  xc  where  the  local  parametrization  (5.4)  is  valid  Then,  if  u  is 
the  solution  of  (5.7b)  satisfying  the  initial  condition  u(tc)=0,  it  follows  that 


(7.2) 


x(tk)  =  0(u(tk))  k=1 ,2 . N,  x(tc)  =  0(u(tc))  =  xc. 
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For  the  approximate  solution  of  (5.7b)  we  consider  a  multistep  method  of 
the  general  form 

(7.3)  uk  =  £  ajUk.j  +  £  pjU’k-i ,  uVi  =  Qi  Ti(0(uk.j)) 

i-1  i-0 

which  is  assumed  to  be  consistent;  that  is, 

(7.4)  i>  =  1, 

i-i 

as  well  as  convergent  of  order  v  >1 .  Assume  that  for  some  k,  1  <  k  <  N,  the 
approximations 

(7.5)  Uj  =  u(tj)  ,  Xj  s  <D(Uj )  ,  j=0,1  ,...k-1 , 

have  already  been  computed.  Then  the  next  point  uk=  u(tk)  along  the  solution 
of  (5.7b)  is  specified  by  (7.3)  and,  by  solving  the  equation  (6.1)  with  u  =  uk,  we 
obtain  xk  =  <X>(uk )  as  the  next  approximation  of  the  solution  of  (4.7). 

As  shown  in  section  5,  the  vectors  u'k_j  of  (7.3)  are  the  solutions  of  a  linear 
system  of  the  form  (5.8).  We  consider  first  the  case  p0  =  0  of  explicit  integration. 
Since  the  local  coordinate  map  (5.4)  is  defined  by  the  matrices  we  have, 
by  (6.1), 

(7.6)  u  =  Q1*(x*xc) 

In  other  words,  we  may  assume  that  the  computed  approximations  (7.4)  satisfy 

(7.7)  uj  =  Q1*(xj-xc),  j  =  1,2 . k-1 

Together  with  (7.6)  this  implies  that  when  uk  is  determined  by  the  multistep 
method  (7.3)  then,  instead  of  using  (6.1)  with  u  =  uk,  we  may  obtain  xk=  0(uk) 
as  a  solution  of  the  nonlinear  system 
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(7.8) 


Fi(x)-b 


A 


=  0 


V 


[GCjXk-i  +  hPm(xk.j)] 


J 


This  leads  to  the  following  algorithm  for  computing  from  the  last  p  known 
approximations 


(7.9)  xN  =  x(tk.j),  2k.j  =  z(tM  ),  wk.j  =  ri(xk.j),j=1,2 . p 

a  next  approximate  point  on  the  desired  solution  of  (4.7): 

Algorithm  !.  (1 )  Evaluate  ak  =  Q*i(^  aiXk.j  +  h^  PjWk.j) 

i-i  j.i 

(2)  Compute  the  solution  xk  of  the  nonlinear  system 

'  FiM-b  N|  =  0 
,  Qix-ak  , 


(3)  With  x  =  x  k  solve  the  linear  system  (5.8)  to  obtain 
wk  =  Q1T]1  and  z  k  =  z. 

This  algorithm  may  be  applied  as  long  as  the  computed  points  do  not  leave 
the  region  of  validity  of  the  local  parametrization  at  xc  defined  by  Q1  and  Q2.  To 
determine  when  the  points  leave  this  region  is  a  very  delicate  problem  and 
requires  special  attention.  In  our  numerical  implementation  we  decide  to 
change  the  local  parametrization  when  the  estimated  condition  number  of  the 
Jacobian  of  the  system  in  step  (2)  becomes  too  large  or  when  the  number  of 
Newton  steps  required  for  solving  this  system  exceeds  a  certain  bound  .  Of 
course,  if  the  multistep  method  (7.3)  has  order  v  then  one  should  solve  the 
systems  in  step  (2)  and  (3)  with  an  error  not  exceeding  0(hv>+1). 


A  further  delicate  problem  is  the  initialization  of  the  overall  process.  One 
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approach  is  to  obtain  u1fu2 . up  by  means  of  some  known  start-up  process  for 

the  multistep  method  (7.3)  applied  to  the  ODE  (5.7b)  and  then  to  set 

Xj  =  0(Uj ) ,  wj  =  ri(xi)  ,i  =  1,2 . p 

by  solving  the  corresponding  systems  (6.1)  and  (5.8)  .  Numerical  experience 
indicates  that  no  reinitialization  is  necessary  when  passing  from  one  local 
parametrization  to  another. 

In  the  case  (30  *  0  the  new  point  is  implicitely  defined  by  both  (7.3)  and 
(5.8);  that  is,  the  systems  in  step  (2)  and  (3)  of  algorithm  1  can  no  longer  be 
solved  separately.  Thus,  in  this  case  the  algorithm  has  to  be  modified  as 
follows: 


Algorithm  2. 


(1)  Evaluate  ai< 


=q;<£ 


ajXk-i 


(2)  Obtain  xk  =  x,zk  =  z,wk  =  was  the  numerical  solution  of 
the  nonlinear  system 


F-i(x)  -  b 

Qi(x-h(30w)-ak 
DF^xJQtW  +  DFi(x)Q2v 
v  A(x)QiW  +  A(x)Q2v  +  B(x)z  +  G(x) 


=  0 


The  comments  following  algorithm  1  apply  again;  in  particular,  the 
nonlinear  system  in  step  (2)  has  to  be  solved  with  an  error  not  exceeding 
0(hu+1). 


These  results  can  be  readily  carried  over  to  second  order  systems  of  the 
form  (4.17).  For  this  we  assume  now  that  the  hypotheses  of  Theorem  6  are 
satisfied  and  hence  that  the  existence  of  a  solution  x:  J  ->  Ex ,  z:J  ->  Ez  of 
(4.17)  is  guaranteed  which  satisfies  the  initial  conditions  x(0)  =  x0,  x'(0)  =  y0, 
x"(0)=r|(x0,y0),  z(0)  =  C(x0,y0)  for  a  given  point  (x0,y0)  e  Ex  x  (Eyn  kerOF^)) 
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As  above,  let  xc  =  x(tc)  e  Mb  b=F-[(x0),  be  a  current  point  on  this  solution 
where  we  construct  a  local  parametrization  (5.4)  defined  by  the  linear  maps 
QvQ2.  In  section  5  we  saw  that  in  some  neighborhood  of  xc,  where  this 
parametrization  is  valid  ,  the  functions 

(7.10)  u  =  Q1*(x-xc),  v  =  Q/x*  . 

satisfy  the  ODE  (5.14). 

In  the  second  order  case,  the  multistep  method  (7.3)  now  reads 

(7.11a)  uk  =  jr<xiuk.i  + 

i-1  i=0  i-1  i-0 


PiV'k-i .  Vk  =  ^  aiVk.j  +  ^ 


PiV'K-j 


where 

(7. 1 1  b)  v'j  =  Q/n(  <D(u-\  D<P(u()vf  ). 

As  we  saw  at  the  end  of  section  5,  Vj'  is  obtained  by  solving  the  linear  system 
(5.15). 

We  begin  again  with  the  case  p0  =  0  and  use  the  same  notation  and 
assumptions  as  before.  In  particular,  let 

(7.12)  xj  =  x(tj),  yj  =  x’(tj),  Wj  =  ri(Xi,  y, ),  z\  =  z(tj) ,  i  =  1,2 . k-1 

be  known  approximations.  Then  we  obtain  the  following  the  explicit-integration 
algorithm  for  computing  from  the  last  p  points  of  this  sequence  the  next 
approximate  point: 

Algorithm-^ 

(1)  Evaluate  a*  =  Q*i(J],  ajXk-i  +  h 

i-1  i-1  i-1  i-1 


Piyk-i).  S'k 


-out 


cxiyk-i  +  h£  pjWk.i) 
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(2)  Solve  the  nonlinear  system 

'  F,(x)  -  b  ' 
Qix*ak  =Q 
Qiy  -  a’k 
v  DF^xJy  , 

to  obtain  xk=x,  yk=y. 


(3)  With  x=xk,  y=yk  obtain  wk  =  Q1rj1  and  zk  =  z  as  solutions  of  the  linear 
system  (5.15). 

When  p0  *  0  then,  analogous  to  algorithm  2,  we  can  implement  the 
multistep  method  (7.1  la/b)  for  the  numerical  solution  of  (4.17)  as  follows: 

Algorithm  4. 

(1)  Evaluate  ak,a’k  as  in  step  (1)  of  algorithm  3. 

(2)  Solve  the  nonlinear  system 

'  F,(x)  -  b 

Qi(x-pohy)  -  ak 

Qi(y-pohw)- a’k  _Q 

DFi(x)y 

DF^xJQ^i+DFiMCfeTte  +  D2F1(x)(y,y) 

,  A(x,y)Q1TH+A(xIy)Q2Ti2+B(x,y)z  +  G(x,y)  , 

to  obtain  xk=x,  yk=y,  wk=  Qirj-j ,  zk=z. 

Obviously  the  comments  made  immediately  following  algorithm  1  apply 
again  to  both  algorithm  3  and  4. 
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8.  Application  to  the  Euler-Laaranae  Equations. 

A  classical  example  for  second  order  DAEs  of  the  form  (4.1 7)  are,  of 
course,  the  equations  arising  in  constrained  multibody  dynamics.  These 
equations  have  the  general  form 

(8.1)  mt)  =0, 

M(q,t)  q"  +  DqvF(q,t)Tz  +  K(q,q\t)  =  0, 


and  are  subsumed  under  the  special  case  of  (4.17)  when  both  A  and  B 
depend  only  on  x.  In  fact,  we  need  to  set  only  x=(q,t)  and  define 


(8.2) 


F,(x)  =  %t),  A(i)=(“M  ° 


r 


B(x)  = 


Dq'Ffx.t) 
0 


T  A 


,  G(x.x')  =  (  K(^q,>t)  ) 


Suppose  that  in  (8.1)  the  mappings'?:  R^xR1-*  Rr,  K^R"*1)2*1^1-*  Rn_1, 
and  MiR^xR1-*  LtR^.R0’1),  r  <  n-1 ,  satisfy,  for  all  (q,t)  under  consideration, 
the  smoothness  conditions  corresponding  to  (4.12a).  Evidently,  (4.12b)  is  here 
equivalent  with 


(8.3)  rankDq'?(q,t)  =  r,  aTM(q,t)a  *  0,  V  a  e  ker  Dq'?(q,t); 


that  is,  with  the  assumptions  that  the  algebraic  constraints  are  independent 
and  that  the  mass  matrix  M  is  definite  on  the  nullspace  of  Dq1?.  Thus,  Theorem 
6  applies  to  (8.1).  Note  also  that  now  the  defining  relations  (4.19)  for  the 
mappings  ri^R""1)2  — >  Rn_1  ,  (Rn_1)2  -»  Rr,  have  the  form 


(8.4) 

f  Ti(q.q’.t)  "i 
l  C(q.q’.t)  J 


M(q,t) 

,  Dq'?(q,t) 


DqmoT  T 

'  (  K(q,qM) 

o  , 

,  D^r(q.t)(q',q')+2D?,4'(q,t)q'+D?,T(u,t)  , 
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We  reduce  {8.1 )  to  a  system  of  ODEs  by  using  a  local  parametrizations  of 
the  type  considered  in  section  6.  However  in  doing  so  care  should  be  taken 
not  to  modify  the  "special"  variable  t .  By  (8.2)  we  have 

(8.5)  DF1(x)  =  (Dq^(q,t),Dt'F(qtt)) 

and  because  of  (8.3)  there  exist  always  r  linearly  independent  columns  among 
the  first  n-1  columns  of  DF^x).  Hence,  when  we  construct  a  local 
parametrization  by  "coordinate  partitioning"  at  a  current  point  xc  =  (qc,tc)  then 

we  can  always  choose  a  permutation  ji,...,jn-i  of  the  set  (1,2 . n-1}  such  that 

the  matrices 


(8-6)  Q.  =  (et. . ew„  en)  =  (£  °  ).  02-  (ej, . 9j)  -  (^) 

satisfy  (5.3)  and  (6.2)  and  hence  define  the  desired  parametrization.  Evidently, 
in  (8.6)  the  (n-l)xr  matrix  P1  and  the  (n-l)x(n-l-r)  matrix  P2  consist  of  the  basis 

vectors  of  R"*1  with  indices  j1(...,jr  and  jr+1 . jn_.,,  respectively.  In  other  words, 

we  construct  the  parametrization  (8.6)  by  working  only  with  the  rx(n-1 )  matrix 
Dq^(qc,tc)  rather  than  with  the  rxn  matrix  DF^xJ. 

Analogously  we  can  implement  a  "tangent  space  parametrization  "  by 
performing  a  QR-factorization  of  Dq'F(qctc)*  rather  than  of  DF^xJ  .  More 
specifically,  if 


(8.7)  Dq'P(qe,tc)’  -  P  (  *  ) 

where  P  is  an  (n-l)-dimensional  unitary  matrix  while  R  is  again  rxr  upper- 
triangular,  then,  corresponding  to  (6.6),  the  matrices 


(8.8)  Qi=(P0)er+1 . en-i,en)  = 


satisfy  (5.3)  and  (6.2)  and  hence  define  the  desired  parametrization.  In  (8.8), 
P2  is  the  (n-l)xr  matrix  consisting  of  the  first  r  columns  of  P,  and  the 
remaining  (n-1-r)  columns  of  P  make  up  the  (n-l)x(n-l-r)  matrix  Pv 
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Now,  the  algorithms  3  and  4  of  section  can  be  applied.  However  since  the 
matrices  Q:  ,Q2  of  (8.6)  or  (8.8)  have  a  special  form  some  simplifications  are 
possible.  Moreover,  it  should  be  taken  into  account  that  in  many  applications 
involving  the  equations  (8.1)  the  accelerations 

(8.9)  x"  =  (o")“'n(X'X') 

are  explicitly  required. 

Suppose  again  that  the  approximate  points 

Xj  =  (Qj.  tj),  y j  =  (q'j.1 ),  Wj  =  (q"j,0),  Z; ,  tj  =  ih,  i=0,1 . k-1 

are  already  known  and  that  xc  =  (qc,tc)  is  the  current  point  for  the  local 
parametrization.  Note  that  the  multistep  method  (7.1  la/b)  was  assumed  to  be 
consistent  so  that,  besides  (7.4),  the  coefficients  must  satisfy 

j-1  j-0 


Hence,  in  step  (2)  of  algorithm  3  and  4,  the  last  equation  of  Q/x  =  reduces 
to  t  =  kh  while  the  last  equation  of  Q,*y  *  a'k  is  trivially  satisfied.  Thus,  the 
explicit  algorithm  3  can  be  simplified  as  follows: 

Algorithm  5. 

ctjqk-i  +  hj)  piq'k-i),  s'k  =  Pi 

i-1  i-1  n1  i-1 

(2)  Obtain  qk=q,  q'k  =q'  as  solutions  of  the  nonlinear  system 


(1 )  Evaluate  ak 


otjq'k-i  +  hX  pfl’k-j) 
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mu 

piq  -  ak  =0 
P‘iq'  -  a'k 

,  DqmtkXf +  Dt^(q.tk)  , 

(3)  With  q=qk,  q'=q’k  solve  the  linear  system  (8.4)  to  obtain  q"k  =  qM, 
and  \  =  z 

For  the  implicit  algorithm,  it  is  useful  to  introduce  the  abbreviation 
Y(q.q'.t)  =  •DqqvF(q,t)(q',q')  -  2D^r(q,t)q'  -  D?t>P(q,t). 


Then,  with  the  acceleration  (8.9),  the  last  equations  in  the  nonlinear  system  of 
step  2  of  algorithm  4  can  be  written  in  the  form 


Dqm  tk)q"  =  Y(q.  q’.  tk) 

M(q,  tk)q"  +  Dq»F(q,  tk)Tz  +  K(q,  q\  tk)  =  0 

Evidently,  the  first  of  these  equations  can  also  be  obtained  formally  by 
differentiating  the  first  equation  of  (8.1)  twice. 

With  these  observations  ,  algorithms  4  reduces  to  the  following  algorithm: 

Algorithm  6. 

(1)  Evaluate  ak,a'k  as  in  step  (1)  of  algorithm  5. 

(2)  Obtain  qk=q,  q'k=q',  q"k=  q",  zk=z  as  solution  of  the  nonlinear  system 

'  m  tk) 

Pi(q-pohq’)  -  ak 
Pi(q'-p0hq")  -  a'k  _q 

Dq4/(q,  tk)q’  +  D,^(q,  tk) 

DqxF(q,  tk)q"  +  y(q,q\tk) 

,M(q,tk)q"+Dq'P(q,  tk)Tz  +  K(q,q',tk) , 


31 


The  numerical  solution  of  the  nonlinear  system  (9.23)  represents  a 
significant  computational  challenge.  However,  in  a  separate  paper  it  shall  be 
shown  that  when  proper  care  is  taken  to  exploit  its  special  structure,  then  the 
computational  complexity  of  this  problem  can  be  reduced  greatly. 


10.  A  numerical  example 

A  frequently  encountered  basic  mechanisms  is  the  four-bar  system  shown 
schematically  in  Figure  1 .  To  model  this  system,  each  of  the  four  links  and  the 
ground  are  considered  as  a  body  and  marked  by  the  integers  1 , 2,  3,  4.  A 
frame  is  fixed  for  each  body  with  the  origin  at  its  centroid.  For  each  body  the 
coordinates  of  the  centroid  are  given  in  Table  1.  This  identifies  also  the  initial 
position  of  the  system  as  it  is  sketched  in  Figure  1 .  Table  2  lists  the  mass  and 
moment  of  inertia  for  each  body.  All  the  joints  of  the  mechanism  are  revolute 
joints,  they  are  marked  by  the  letters  A,  B,  C,  D.  Table  3  identifies  the  two 
bodies  connected  at  each  joint  and  gives  the  coordinates  of  the  common  point 
in  the  coordinate  frame  of  each  body.  .  The  mechanism  is  subject  to  no 
applied  force  other  than  gravity.  The  initial  velocity  of  tuch  generalized 
coordinate  is  assumed  to  be  equal  to  zero. 


Body  # 

1 

2 

3 

4 

X 

0.0 

2.0 

6.0 

12.0 

y 

0.0 

2.0 

4.0 

4.0 

0 

0.0 

7t  /  4 

0.0 

71/4 

Table  1  :  Body  centroids 


Body  #1  2  3  4 


Mass  (kg)  1.0  200.0  300.0  100.0 

Inertia  (kgm)2  1.0  200.0  300.0  100.0 


Table  2:  Masses  and  Moments  of  Inertia 
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Joint 

A 

B 

C 

D 

Bodyi 

1 

2 

3 

4 

Xj(P) 

0.0 

2.828 

4.0 

2.0 

yj(P) 

0.0 

0.0 

0.0 

0.0 

Body  j 

2 

3 

4 

1 

Xj(P) 

-2.828 

-4.0 

-2.0 

12.0 

Vj(P) 

0.0 

0.0 

0.0 

0.0 

Table  3:  Joint  data  at  the  common  point  P 


The  constrained  equations  of  motion  for  this  system  have  been  solved 
using  the  algorithms  discussed  in  Section  9.  At  the  same  time,  as  a  check  the 
model  was  also  analyzed  by  means  of  DADS,  a  general  purpose  code  for 
mechanical  design  [2]. 

The  numerical  experiments  were  carried  out  with  the  implicit  as  well  as  the 
explicit  algorithm.  The  second  order  BDF-formula  was  chosen  as  the  implicit 
method  and  the  second  order  PECE-  (Adams  -Bashforth  -  Moulton)-formula  as 
the  explicit  method.  In  each  case,  fixed-stepsizes  h=0.01 ,  and  h=0.0025  were 
applied.  The  nonlinear  systems  were  solved  by  a  chord  Newton  method.  As 
indicated  earlier,  the  local  parametrization  was  changed  whenever  the 
estimated  condition  number  of  the  linear  systems  exceeded  a  specific  bound 
or  when  the  nonlinear  solver  failed  to  converge  in  a  given  number  of  steps.  In 
Figures  2-4  we  show  the  position,  velocity,  and  acceleration  of  the  x- 
coordinate  of  the  centroid  of  body  2  for  the  tangent  space  parametrization. 


0.  0.325  0.65  0.975  1.3  1.625  1.95  2.275  2.6  2.925  3.25 

x  E  0 

TIME  (SEC) 


Figure  2:  Position  of  the  x-coordinate  of  the  centroid  of  body  2. 

x  E  1 


o.a 

0.575 

0.35 

0.125 

-0.1 

-0.325 

-0.55 

-0.775 

-1. 

0.  0.325  0.65  0.975  1.3  1.625  1.95  2.275  2.6  2.925  3.25 

x  E  0 


TIME  [SEC) 


Figure  3.  Velocity  of  the  x-coordinate  of  the  centroid  of  body  2. 
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0.  0.325  0  .  j.975  1.3  1.625  1.95  2.275  2.6  2.925  3.25 

x  E  0 

TIME  (SEC) 


Figure  4:  Acceleration  of  the  x-coordinate  of  the  centroid  of  body  1 . 


The  numerical  solutions  obtained  with  the  implicit  method  are  accurate  to 
four  decimal  digits  throughout  the  time  interval  [0,3]  and  hence,  in  the  figures, 
they  cannot  be  distinguished  from  the  solutions  given  by  DADS.  On  the  other 
hand,  the  solutions  obtained  with  the  explicit  method  are  drifting  away  for  t  > 
2.5  sec.  (see  Fig  2,  3,  4)  when  the  the  stepsize  is  increased. 

The  computations  were  also  carried  out  with  the  local  parametrization, 
induced  by  "coordinate  partitioning  "  and  they  agree  to  three  significant  digits 
with  the  result  for  the  "tangent  space"  parametrization.  Table  4  provides  some 
statistical  information  about  the  total  number  of  function  and  Jacobian 
evaluations  as  well  as  about  the  number  of  re-parametrizations  for  the  two 
choices  of  local  parametrizations  when  the  step-size  h  =  0.01  was  used.  In  the 
table,  TG  and  GCP  stands  for  the  tangent  space  and  coordinate  partitioning 
parametrization,  respectively. 
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Explicit  Method 

Implicit  Method 

TG 

GCP 

TG 

GCP 

Function  Eval. 

989 

1021 

1456 

1520 

Jacobian  Eval. 

56 

65 

134 

153 

Re-Parametr. 

6 

6 

15 

19 

Table  4 : 

Performance  Statistics 
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